arXiv: 1506.01137v2 [q-bio.PE] 18Jun2015 


Stochastic dynamics and logistic population growth 


Viceng Mendez/ Michael Assaf/ Daniel Campos/ and Werner Horsthemke^ 

^ Grup de Fisica Estadistica. Departament de Fisica. Facultat de Ciencies. Edifici 
Cc. Universitat Autonoma de Barcelona, 08193 Bellaterra (Barcelona) Spain 
^Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel 
^ Grup de Fisica Estadistica. Departament de Fisica. Facultat de Ciencies. Edifici 
Cc. Universitat Autdnoma de Barcelona, 08193 Bellaterra (Barcelona) Spain 
'^Department of Chemistry, Southern Methodist University, Dallas, Texas 75275-0314, USA 

(Dated: June 19, 2015) 

The Verhulst model is probably the best known macroscopic rate equation in population ecology. 
It depends on two parameters, the intrinsic growth rate and the carrying capacity. These parameters 
can be estimated for different populations and are related to the reproductive fitness and the com¬ 
petition for limited resources, respectively. We investigate analytically and numerically the simplest 
possible microscopic scenarios that give rise to the logistic equation in the deterministic mean-field 
limit. We provide a definition of the two parameters of the Verhulst equation in terms of microscopic 
parameters. In addition, we derive the conditions for extinction or persistence of the population 
by employing either the “momentum-space” spectral theory or the “real-space” Wentzel- Kramers- 
Brillouin (WKB) approximation to determine the probability distribution function and the mean 
time to extinction of the population. Our analytical results agree well with numerical simulations. 

PACS numbers: 05.40.-a,87.23.Cc,87.10.Mn 


I. INTRODUCTION 


Quantitative models of population dynamics have at¬ 
tracted an enormous interest from biology to mathemat¬ 
ics and physics [H-Q. In the deterministic limit, these 
models coincide with macroscopic rate equations based 
on phenomenological laws. The simplest one corresponds 
to Malthus law, where the per capita rate of change in 
the number of individuals is constant, resulting in a linear 
growth rate for the population, dn/dt = rn. The popu¬ 
lation grows exponentially, n{t) = noexp(rt), where r is 
the intrinsic growth rate and rig is the initial population. 
Unlimited exponential growth is patently unrealistic, and 
factors that regulate growth must be taken into account. 
The most famous extension of the exponential growth 
model is the Verhulst model, also known as the logistic 
model, where the per capita rate of change decreases lin¬ 
early with the population size. The population’s growth 
rate, dn/dt = rn(l — n/K), is now a quadratic function 
of the population size, where K is known as the car¬ 
rying capacity. This equation was derived initially by 
P. Verhulst in 1845 [3, 0 and was rediscovered later by 
R. Pearl in 1920 Other models, like the Gompertz 
growth, dn/dt = anln{K/n), exhibit many of the same 
properties, but the logistic equation is arguably the best- 
known and most widely applied rate equation for popu¬ 
lation growth and population invasion [l|, l7| . 

These models are deterministic and ignore fluctua¬ 
tions. Real populations evolve in a stochastic manner, ex¬ 
periencing intrinsic noise (or internal fluctuations) caused 
by the discreteness of individuals and the stochastic na¬ 
ture of their interactions, see, e.g., [Ml- When the 
typical size of the population is large, fluctuations in 
the observed number of individuals are typically small 
in the absence of external or environmental noise. The 


dynamics of the population then can be described by a 
deterministic mean-field rate equation. In the case of the 
logistic equation, the population evolves from an initial 
condition to a stable stationary state, where the popu¬ 
lation size equals the carrying capacity and persists for¬ 
ever. However, if the typical population size is not large, 
internal fluctuations can lead to the extinction of the 
population [l^. The effects of internal fluctuations have 
been studied in predator-prey models EH, epidemic 
models fl7l - E3 . cell biology and ecological systems 


In particular, extinction of a stochastic population 
11 , [IS , which is a crucial concern for population bi¬ 


ology w and epidemiology jz 


I, has also attracted 

scrutiny in cell biochemistry [30[ and in physics 


To describe the intrinsic noise of populations, we adopt 
individual-based models, also called stochastic single 
patch models [H, . An individual-based formulation 

provides several advantages. It is often easier to define 
an ecological system in terms of the events that govern 
the dynamics of the system at the level of individuals. 
Population-level models, such as the Verhulst equation, 
can then be derived analytically as the mean-field approx¬ 
imation, instead of simply be postulated phenomenolog¬ 
ically. In this way, individual-based models provide a 
microscopic basis for the usual ecological rate equations, 
and the range validity of the latter can be established 
by comparing its predictions with those of the former. 
Individual-based models capture the fact that popula¬ 
tions consist of discrete individuals undergoing random 
events corresponding to birth (reproduction of the popu¬ 
lation), competition (between individuals for limited re¬ 
sources) and death (natural decay of individuals). It is 
well known that different types of of individual-based 
schemes are described by the same Verhulst equation in 
the deterministic limit. Since extinction is ultimately 
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caused by the stochastic nature of the interactions be¬ 
tween individuals, it is critical to analyze how the details 
of the individual processes affect the ultimate fate of the 
population or the time to extinction. We explore a vari¬ 
ety of stochastic interactions between individuals, all of 
which give rise to the logistic equation in the mean-field 
limit. We find different dynamical behaviors, such as per¬ 
sistence or extinction, of a population that experiences 
birth, death, and competition processes. Extinction is 
due to rare fluctuations, and the mean extinction time 
(MET) of the population strongly depends on the mi¬ 
croscopic details of the processes, such as the number 
of “newborn” individuals or the number of individuals 
removed due to exclusive competition. We obtain ana¬ 
lytical solutions for the probability distribution function 
(PDF) of individuals, if the population persists, and for 
the MET, if the population becomes extinct. Our ana¬ 
lytical results are compared with numerical simulations, 
performed using the first reaction method [s^. We con¬ 
sider the birth-and-death and birth-competition-death 
cases separately, making use of the “momentum space” 
spectral theory [^,[^ and the “real-space” WKB theory 
[3§ - l40l |. respectively. 

II. MASTER AND MEAN FIELD EQUATIONS 
FOR GENERAL BIRTH-COMPETITION-DEATH 
PROCESSES 

We investigate individual-based models of populations 
in which the following birth, competition, and death pro¬ 
cesses occur, 


bX\{a + b)X, (2.1a) 

cXA(c-(i)X, (2.1b) 

X^0, (2.1c) 

where a, 6, c, and d are positive integers, and d < c. Such 
processes occur also in chemically reacting systems, and 
it is convenient to adopt the language of chemical kinet¬ 
ics to make a connection with the literature of stochas¬ 
tic chemical models. Therefore, we will often refer to 
the processes of (12.11) as “reactions.” If d = c, the last 
two reactions are death reactions, due to competition 
between c individuals {cX 0) or due to natural de¬ 
cay {X A 0). We make the standard assumption that 
the reaction scheme (12.11) defines a Markovian birth-and- 
death process, see, e.g., [IE IHi EH) Ell j and employ the 
Master equation, also known as the forward Kolmogorov 
equation, to describe the temporal evolution of P{n,t), 
the probability of having n individuals at time t, 

~ ^ ^~ — r,t) — W (n, r)P{n, t)]. 

r 

( 2 . 2 ) 

Here, IE(n, r) are the transition rates between the states 
with n and n -I- r individuals, and r = {ri,r 2 ,r 3 } = 


{a, —d, —1} are the transition increments. Equation (12.2p 
can generally only be solved in the stationary limit, 
dP{n,t)/dt = 0 and only for the special case that 
a = d = 1, i.e., only single-step processes occur in the 
population. Then the condition of detailed balance holds, 
which significantly simplifies the theoretical analysis, and 
exact analytical expressions can in principle be obtained 
for the stationary PDF or the MET. For recent reviews, 
see for example EE Ej) E3 ■ We emphasize that we study 
the general generic case of arbitrary a and d to elucidate 
how the microscopic details affect the PDF or the MET. 
Calculating the stationary PDF or the MET is highly 
nontrivial for multi-step reactions, and this case has only 
recently began to be addressed. 

The transition rates corresponding to each reaction, 
W(n,r), are obtained from the reaction kinetics and 
for (12.ip read: 

W (n,—1) = jn. (2.3c) 

Substituting (12.31) into (j2.2l) . we find 


dP{n,t) X (n —a)! 


dt b\{n-a-b)\ 

^ {n + dy. 


P(n — a, t) 


c! {n + d — c)\ 
A n! 


P(n + d,t) + 7 (n -|- l)P{n + 1, t) 


fj, n\ 


bl (n — by cl (n — c)! 


■ yn 


P(n,t), (2.4) 


where it is understood that P{n <0,t) = 0. The proba¬ 
bility generating function [s^ is defined as 


G'(p, i) = X! (2-5) 

n—0 


where p is an auxiliary variable, which is conjugate to 
the number of particles EE- Once G{p,t) is known, the 
PDF is given by the Taylor coefficients 


P{n,t) = — 

n! 


5"G(p, t) 




( 2 . 6 ) 


- p—Q 


Normalization of P{n,t) implies that G{p = l,t) = 1. 
Multiplying (12.41) by p", summing over n, and renaming 
the index of summation, we find 




dt 


bl 


n—0 

CO 


(n — by.' 


{n-cy 

CO 

+ lY.{p^-^-pnnPM. (2.7) 


n=0 
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Taking into account the property 

^ n(n - 1) • • • (n - /c + l)p”P(n, t) 

^ n^O 

oo I 

= P-8) 

n=0 ^ ^ 

in (lO) . we finally obtain the evolution equation for 


DGjP.t) ^ A , _ ^ 

at ’ dpb 




Equation (1^ is exact and equivalent to the Master 
equation (12.41) . If only one individual reactant is present 
in all the reactions, i.e., b = c = 1, then (1^ is first 
order in p and can be solved exactly using the method of 
characteristics. 

Macroscopic equations, i.e., equations for the expected 
or average values, can be obtained easily from (12.411 . Mul¬ 
tiplying (12.41) by , summing up over n, and renaming 
the index of summation, we find 


dt 




n—0 


n=0 ^ ^ 

oo I 


n=0 


+ -f^[in-l)'^-n'^]nPin,t). ( 2 . 10 ) 


n—0 


The fc-th moment is defined as (n*) = Yl^=on^P{n,t) 
and evolves according to the ordinary differential equa¬ 
tion 


where p = (n) is a macroscopic quantity, the average or 
expected number of individuals in the population. 


III. BIRTH AND DEATH/COMPETITION 
PROCESSES 

We consider the case of two reactions, i.e., 7 = 0: 

6X^(&-ha)X, (3.1a) 

cXA(c-d)X. (3.1b) 

In the first reaction (birth), b individuals have to interact 
with each other to produce a new individuals at a con¬ 
stant rate A. In the second reaction (death by compe¬ 
tition), c individuals interact with each other to remove 
d individuals at a constant rate p. The fact that b, in 
general, can be larger than 1 includes scenarios where a 
single individual cannot generate by itself new individu¬ 
als, which represents a type of Allee effect (i^ . 

Equation (12.121) reduces to the logistic equation if & = 1 
and c = 2. From a kinetic point of view this means 
that an individual does not need to interact to give rise 
to new individuals; the birth reaction takes the form 

X A (a -I- 1) X. The fact c = 2 implies that a linear 
death rate, corresponding to X —>• 0, cannot occur for 

the scheme (ED in this case. The possible death re¬ 

actions, compatible with a mean-field logistic equation, 
are 2 X A X (competition) or 2 X A 0 (annihilation). 
Consequently, the birth-and-death processes that lead to 
logistic macroscopic behavior are 

xA(a-bl)X, (3.2a) 

2XAx, (3.2b) 

and 

xA(a-bl)X, (3.3a) 

2XA0. (3.3b) 


d{n^) 

dt 


6-1 


6 ! 


[(n -I- a)^ — n^] (n — m) 


m—0 

c-1 


+ 3 ( [A - d)'" - n'^] n 


m—O 


+ 7 ( [(n - 1)'= - «'=] n). 


( 2 . 11 ) 


The logistic equation for these two reaction schemes reads 



where 


r = aXa,nd N = 2aX/pd (3-5) 


Equation (12.111) is not closed, and one must deal with 
a hierarchy of coupled differential equations for k = 
1,2,3,.... In order to truncate this set and to obtain 
closed equations, we make use of the mean-field approxi¬ 
mation ~ (n)^. which holds if the typical population 
size is large [HI, [s^- For k = 1, the mean-field equation 
reads 


dp 

dt 


= ^p^ 

b\ ^ 


pd 


IP, 


( 2 . 12 ) 


are the intrinsic growth rate and the carrying capacity, 
respectively. These definitions are valuable because they 
allow us to relate the macroscopic parameters r and N, 
which can be measured for different kinds of popula¬ 
tions, to the microscopic parameters that characterize the 
stochastic processes involved in the interaction between 
the individuals of the population. From a macroscopic 
point of view, the logistic equation for population growth 
is specified by two parameters. On the other hand, the 
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schemes (IQ) and (IQ) contain up to four microscopic 
parameters, namely a, d, fi and A. As a result, we have 
two additional free microscopic parameters that can take 
arbitrary positive values compatible with the same mean- 
field logistic equation. Rate equation (13.41) has an unsta¬ 
ble steady state at ps = 0 and a stable steady state at 
P;, = for d = 1 or d = 2. Below we deal separately with 
schemes dm and dm and apply the momentum-space 
(p-space) spectral theory to find the stationary PDF in 
the case of population survival or the MET in the case 
of population extinction. An important advantage of the 
p-space representation stems from the fact that the evo¬ 
lution equation for the generating function G'(p, t) is ex¬ 
actly equivalent to the original master equation. There¬ 
fore the p-space approach is especially valuable for an 
exact analysis. 


A. Case I: X ^ (a-t 1)X, 2X A X 


In this case we expect the population to evolve to a 
nontrivial steady state and not to become extinct. The 
equation for the probability generating function, (1^ . 
becomes 


dGip,t) ^ , p, 2^3^G' 


(3.6) 


If initially at t = 0 the system consists of ng individuals, 
then P{n, 0) = Sn,no, where S is the Kronecker delta, and 
from (|2.5I) we find G(p, t = 0) = . The boundary con¬ 

ditions (BCs) are “self-generated”. Indeed, the equality 
G(p = I,t) = I holds at all times, due to the conserva¬ 
tion of probability. Equation (13.61) has a singular point 
at p = 0. Since G(p,t) must be an analytic function at 
p = 0 for all times, we require that G(p = 0,t) = 0. This 
condition stems from the fact that G(p = 0,t) = Po{t), 
and since the population cannot go extinct, the probabil¬ 
ity of extinction vanishes at all times. We are interested 
in the steady state. Then (13.61) turns into 

|(l-p)G" + A(p“-I)G' =0, (3.7) 

which must be solved with the BCs Gs(l) = I and 
Gs(0) = 0. The exact analytical solution reads 

G.(P) = ( 3 ^ 8 ) 

Jg exp[N(f>{s)/a\ds 

where 


(j){s) = — ln(I — s) — 


I — s 


ds = ^ —, (3.9) 


and N = 2aA/p. In the special case where a = 1, the 
exact solution for the generating function can be easily 
obtained from (13.81) and dm, 


GM 


exp{Np) — I 
exp(A^) - 1 ■ 


(3.10) 


Expanding exp(7Vp) around p = 0, we find that for large 
N the stationary PDF follows the Poisson distribution, 


Ps{n) 


exp(-iV) 


(3.11) 


where we have approximated exp(A^) — 1 ~ exp(A^) in the 
denominator. We have performed numerical simulations 
and compared them with (13.111) . Figure [T] shows that 
the agreement becomes better as the typical number of 
individuals N is increased. 



n n 


FIG. 1: Stationary PDF for X A 2 X, 2 X A X. In panel a) 
N = 100, in b) A = 25, in c) A = 200 and in d) A = 50. 
Simulation results (symbols) are based on 3000 realizations of 
the stochastic process up to time 10®. 


If o = 2, the generating function is given by 

erfi (if) - erfi (Viv) ’ 


(3.12) 


where erfi((r) = exp{t^)dt. The PDF can be ob¬ 

tained by substituting (13.121) into (I2.6|l . A comparison 
between the analytical PDF and numerical simulations 
is shown in Fig. [2l and excellent agreement is observed. 

Finally, we can also obtain the mean number of in¬ 
dividuals in the stationary state and its dependence on 
N = 2aA/p by using the definition of G from (j2.5ll . Dif¬ 
ferentiating (13.8p with respect to p and using (13.9p . we 
find 


(n) = G'(l) = , ■ 

fo expl^(/)(s)Jds 

Furthermore, the variance of n satisfies (n^) 
G"(l) -b G'(l) - G'(l)2, and we find 


(3.13) 

(n)2 = 


(n^) - (n)^ = (n}(l + N) - (n)^. (3.14) 
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FIG. 2: Stationary PDF for X ^ 3X, 2X A X. In panel a) 
N = 200, in b) iV = 50, in c) = 400 and in d) iV = 100. 
Simulation (symbols) results are based on 3000 realizations of 
the stochastic process up to time 10®. 


This allows us to determine the coefficient of variation, 
c„, defined as the ratio of the standard deviation to the 
mean, which measures the variability in relation to the 
mean of the population, 


for the mean, the variance, and the coefficient of variation 
if N is large. In that case, the integral in the denomina¬ 
tor of (2.9) can be evaluated by integration by parts for 
Laplace integrals, and we find 


exp 


N 




ds- — \ exp 


N 


Al) 


As a result, the mean value reads 
(n) N ( 1 + a exp 


--m 


. (3.16) 


(3.17) 


and the coefficient of variation is given by 


1 / a 


-A(‘) 


(3.18) 


B. Case II: X A (a-f 1)X, 2X A 0 

In this case, the initial number of individuals and the 
parameter a play a crucial role in determining the ul¬ 
timate fate of the population. Since the death process 
involves two individuals, population extinction is guar¬ 
anteed, regardless of the initial number of individuals, if 
the number of newborn individuals a is odd, i.e., a -I- 1 
is even. In contrast, if a is even, i.e., a -I- 1 is odd, the 
population becomes extinct only if hq is even. 


1. a is even and no is odd 


^ V{n^) - ^ I (1 + N) 

(n) V (n) 


(3.15) 


In Fig. [3] we plot the coefficient of variation Cv obtained 
from numerical simulations (circles) and compare it with 
the theoretical result given by (13.151) . 



We begin by considering the case that a is even. Then 
the birth process preserves the even-odd parity of the 
number of particles. As a result, the population becomes 
eventually extinct if the initial number of individuals uq is 
even. If no is odd, the case considered in this section, the 
population evolves to a nontrivial stationary state. The 
equation for the probability generating function, (1231), is 
given by 


dG{p,t) dG p ^ d'^G 

= -1)—+ -(l-p)_. 


(3.19) 


The boundary condition G{p = 1, t) = 1 still applies, but 
the singular point of (13.191) occurs at p = —1 and not 
at p = 0 as in (1^ . Since G'(p, t) must be analytic at 
p = —1 for all times, we require that G(p = —l,t) = 
(—1)”°. This boundary condition stems from the fact 
that G(p = —l,t) is the sum of all even probabilities 
minus the sum of all odd probabilities [s^. The steady 
state has to be solved by integrating the equation 


FIG. 3: Goeflicient of variation c„ versus N for a = 1 and 
a = 2. The log-log plot in the inset shows that Cv decays like 
Simulation results are based on 3000 realizations of 
the stochastic process up to time 10®. We set p = 2 and d — I 
and vary A. 

It is straightforward to obtain asymptotic expressions 


|(1 - p2)G" + Ap(p“ - 1)G; = 0, (3.20) 

with the boundary conditions Gs(l) = 1 and Gs(—1) = 
(_l)"o_ xiie exact solution reads 

rp 


Gs{p) = Gi 


exp [A^ip(s)/a] ds -|- G 2 , 


(3.21) 
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where 


(/7(s) 


-lii(l-s2)_ 



gO+l 
1 — 


ds, 


(3.22) 


In Fig. [5] we plot the coefficient of variation c„ for the 
cases of a = 2 and a = 4. The mean number of individ¬ 
uals in the steady state, in) = G'(l), can be determined 
from (mSl), 


and N = aX/fi. For no odd, we obtain from the boundary 
conditions 


Cl = - - -, (3.23) 

exp [Nipis)/a] ds 

and 

G 2 = 1 - ^- - -. (3.24) 

exp [A^(/ 3 (s)/a] ds 

As expected, the system reaches a nontrivial stationary 
state with 


Gs(p) — 1 + 2 


exp [A^(^(s)/a] ds 
exp [A^i^(s)/a] ds 


(3.25) 


To be specific, we focus on the case a = 2, that is X A- 
3X, 2X A 0. From (13.221) we obtain (p(s) = s^, and 
from (13.251) 


Gs(p) = 



(3.26) 


The PDF is obtained by substituting (I3.26P into ()2.6I) . In 
Fig. |3]we plot the PDF Pain) for different values of N. 


. . V^expiN/2) 

(^) = - : -- , 

erii{-\/2N 12) 


(3.27) 


and the coefficient of variation is given by (I3.15|) with (n) 
given by (13.2711 . The solid curve corresponds to the ana¬ 
lytical results, and the symbols correspond to numerical 
simulations. The inset again shows that Cy scales like 
iV-l/2. 



N 


FIG. 5: Coefficient of variation versus N for a = 2 and a = 4. 
The inset demonstrates that Cv decays like . Simulation 

results are based on 3000 realizations of the stochastic process 
up to time 10®. We set p. — d — 2 and vary A. 



FIG. 4: Stationary PDF for X A 3X, 2X A 0. In panel a) 
N = 200, in b) A = 40, in c) A = 400 and in d) A = 100. 
Simulation (symbols) results are based on 3000 realizations of 
the stochastic process up to time 10®. 


2. a and no are even 

If no is even, the boundary conditions lead to 

Cl J exp[A(/j(s)/a]ds-I-G 2 = 1 , (3.28) 


and 

Cl J exp[A(/j(s)/a]ds-I-C 2 = 1, (3.29) 

so that Cl = 0 and C 2 = 1. Therefore, Cs(p) = 1 
which describes an empty population state, i.e., extinc¬ 
tion, as < —>■ 00 . To calculate the MET, we employ the 
“momenturn-space” spectral method developed recently 
[H, 113, m, 113 . After a short relaxation time C, which 
corresponds to the deterministic relaxation time of the 
system to the stable stationary state, the population typ¬ 
ically settles into a long-lived metastable state, which 
is encoded by the lowest excited eigenmode ijiip) of the 
probability generating function C(p, <) [ 13 ]. Indeed, for 
t ^ tr, we can write 

C(p, t) = Gsip) - tpip) exp . 


(3.30) 
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Here Ei is the lowest nonzero eigenvalue, r = is 

the mean time to extinction, and Gs{p) = 1. Substituting 
(13.3011 into (I3.6L we obtain 

(1 - - l)il}'{p) = -2Eitp{p), (3.31) 

where n = X/p. Since a is even, the function ^p{p) is also 
an even function. It is therefore sufficient to consider 
the interval 0 < p < 1. Since (n) ^ H, we assume that 
H ^ 1 to find the eigenvalue Ei, which we expect to be 
exponentially small in H. We will proceed by matching 
the asymptotic expansion for the function tpip) in the 
bulk region, 0 < p < 1, namely 'ipb, with ipi, the solution 
in the boundary layer, 1 — p <C 1. We will show that the 
function '0(p) is almost constant everywhere within the 
interval p S [0,1), except in a narrow layer close to p = 1. 
In the bulk we can treat Ei as a perturbative parameter. 
To zero order we set Ei = 0, and the even solution of 
(13.3111 is 1. To account for corrections, we write ip{p) = 
1 + dip, where dip satisfies the differential equation 


In the boundary layer, 1 — p <C 1, we disregard the ex¬ 
ponentially small term Eiip in (13.3111 and integrate the 
resulting equation (1 — p^)'ip"{p) -\- 2Hp(p“ — l)^/i'(p) = 0 
to obtain 

/ p 

exp [H(p(s)] ds, (3.38) 

where we have made use of the boundary condition at 
the boundary layer, i.e., tpi{l) = 0. Equation (13.3811 can 
be rewritten as 

ipi{p) = C exp [f2(p(s)] ds — J exp [f2(p(s)] ds 

^ ^ _ JJexpjSV^y (3.39) 

\ /o exp[Hp(s)]ds/ 

Matching the solutions (13.3711 and (13.3911 . we find Ci = 1 
and the MET, 


sp,” + 2np^ —(3.32) 

1 — p^ 1 — p^ 

whose solution, using Eq. (13.2211 . takes the form 
5p}'{p) = Co exp [H(p(p)] 

-2i5.»p|!i^(p)] (3.33) 

To solve for pj{p), we need to specify two boundary con¬ 
ditions. Setting p = 0 in p.31ll . we obtain ^””(0) = 
—2Eip!{0), or equivalently 6pj''{0) = —2Ei — 2Ei6pj{0). 
On the other hand, from (13.3211 and setting p = 0, we 
find the first boundary condition, Sp}"{0) = —2Ei. This 
condition together with 5p)"{Q) = —2Ei — 2Ei5p){Q) leads 
to the second boundary condition, <5^(0) = 0. The first 
boundary condition implies that ()3.33|1 reduces to 

5p}'{p) =-2Eieyip[Vtp{p)] f (3.34) 

Jo i — s 

which can be integrated together with the second bound¬ 
ary condition to yield 

3V.(ri = -2B. r exp [%(.)]* 

Jo do 1 - 

(3.35) 

Since this solution holds in the bulk region 1 — p ^ 
with H ^ 1, we can approximate the inner integral in 
(13.3511 as follows 


1 — 


Jo 

POO 

~ / exp [—r2ip(tt)] du. (3.36) 
Jo 


Therefore, 


rp coo 

p>b{p) ~ 1 — 2Ei / exp [Hp(s)] ds / exp [—Hp(m)] du. 

Jo Jo 

(3.37) 


2 /* ^ 

T = — / exp [Hp(s)] ds / exp [—r2ip(u)] du. (3.40) 
M Jo Jo 

Since H ^ 1, we can further approximate (13.4011 . The 
function (p{s), given by p.22p . can be expressed as 


a/2 






(3.41) 


for even a, and 


(a-l)/2 2j-|-l 

<d(s) =-21n(l + s) + 2 (3.42) 

1=0 ^ 


for odd a. Since in this subsection we consider the case 
of even a, <p(s) is a polynomial of order a with positive 
coefficients. Therefore, the main contribution of the first 
integral in (13.4011 comes from the region around s = 1. 
Employing the Taylor expansion we find 


exp [r2(p(s)] ds ~ 


exp [Hp(l)] 
Hp'(l) 



exp{H [(p(l) -I- p'(l)(s - 


l)]}ds 


(3.43) 


For the second integral in (13.4011 . the main contribution 
comes from the region around m = 0. To leading order, 
p{u) ~ and 


poo poo 

/ exp [—Hp(m)] dn ~ / exp [—du = 

Jo Jo 




2\f^' 

(3.44) 

Substituting these results into (13.4011 , we obtain a general 
result for the MET for H 1 and any even a. 


r = 


y^exp(HX;“£l y) 


paH3/2 


(3.45) 














































FIG. 6: Mean time to extinction t vs N (panel a)) and vs 

a (panel b)) for reaction X ^ (o + 1) X, 2 X A- 0. Solid 
curves are obtained from (I3.45II . while symbols correspond to 
numerical simulations. We set = 2 and vary A. Simulations 
have been performed up to time 10®. 


As an example, for a = 2 we find 


^/^^exp (n) 

2/rn3/2 ’ 


(3.46) 


which coincides with the result in [s^. We have verified 
the result (13.451) by numerical simulations. In the upper 
panel of Fig. IH we plot r versus N for a = 2 and a = 4, 
and in the lower panel we plot r versus a for different 
values of n. In all these comparisons we obtain excellent 
agreement between theory and simulations. 


3. a is odd 


now the second boundary condition to be = 0, 

where we have neglected the term £’i'0(—1), which is ex¬ 
ponentially small. The final solution for the function tp 
in the bulk region is very similar to the even a case, and 
we find 

ipbip) = 1 - 2£'i f exp [II(/3(s)] ds f du. 

Jo J-l f — u 

(3.47) 

In the boundary layer we obtain exactly the same result 
as (I3.38p . By matching both solutions in the common 
region, we obtain 


M Jo 


exp [f2(/j(s)] ds 


1 — 


To proceed, we employ the approximations (13.361) and 


exp [—n(^(it)] 
1 — 


du J exp [—r2(/j(u)] du 

/ oo poo 

exp [—II(/?(u)] du ~ / exp [—IIu^] du = 


vn' 

(3.49) 


As a result, similar to the even a case, we obtain from 
()3.42|1 the general result for any odd a, 


T = 


/iaII3/2 


exp 


/ “ 2 ^ 

-2nin2-l-2n^ 

\ 1=0 



(3.50) 


For 0 = 1, (13.501) yields 


2A^exp [217(1 — In 2)] 
/iI73/2 


(3.51) 


which coincides with the result in [ 3 ^ . 

In Fig. [7] we verify the result (13.501) for the MET. In 
the upper panel we plot r versus A^ for a = 1 and o = 3. 
The mean time to extinction increases as the number of 
individuals increases, as expected. In the lower panel we 
plot T versus a for relatively low values of 17, and the 
agreement between theory and numerical simulations is 
still fair. 


IV. BIRTH-COMPETITION-DEATH 
PROCESSES 


If a is odd, (13.191) has no other singularity and we have 
only one boundary condition, Gs(l) = 1. As a result, 
Gs{p) = 1, and the population becomes extinct, regard¬ 
less of the value of no- To obtain the MET in this case, 
we start again with (13.191) . Since ipip) is no longer even, 
the bulk region now corresponds to p G [—1,1), and the 
boundary layer is located at 1—p <C 1. In the bulk region 
we impose the boundary condition 6ijj{0) = 0, as in the 
case of even a. However, setting p = — 1 in (13.311) . we find 


We add the death reaction X A 0 to the system of 
birth-competition processes (EH). To obtain a logistic 


equation in the mean-field limit, we consider 6=1 and 
c = 2, leading to the reaction scheme 

xA(a-bl)X, (4.1a) 

2XA(2-d)X, (4.1b) 

X A 0. (4.1c) 
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FIG. 7: Mean time to extinction r vs fl (panel a)) and vs 

a (panel b)) for reaction X A- (o + 1) X, 2 X A- 0. Solid 
curves are obtained from (I3.50II . while symbols correspond to 
numerical simulations. We set p = 2 and vary A. Simulations 
have been performed up to time 10®. 


Here, a > 1 and d = 1 for a birth-competition-death 
system and d = 2 for a birth-annihilation-death system. 
It is straightforward to show that this system always goes 
extinct. We are interested in calculating the MET for 
the general case. Although this can also be done via 
the generating function (p-space theory], we will use the 
“real-space” WKB approximation [^. I38l - l40l | . According 
to (|2.3I) . the transition rates are given by 

W(n,a) = Xn, (4.2a) 

Win,-d) = l^ = lnin-l), (4.2b) 

W{n,—l) = "fn. (4.2c) 

Replacing t by t/j and introducing the rescaled popula¬ 
tion number density q = n/N, where N = 1, the 

transition rates can be rewritten as 

W{n,r) = W{Nq,r) = Nwr{q)+Ur{q) + 0{N~^), (4.3) 
where 


Wa{q) = Roq, 


(4.4a) 


w-d{q) = (4.4b) 

w-i{q) = q. (4.4c) 

Here q, Wr{q), and Ur{q) are 0(1), and 

Ua{q) = u-i{q) = 0, (4.5a) 

u-d{q) = -^Roq- (4.5b) 

Further, Rq = A/y is the basic reproductive number. 
Since n = g = 0 is an absorbing state (extinction), we 
have Wr{0) = Ur.(0) = 0 for any r = jg, —d, —1|. For 
1, the WKB theory developed in (^ l38M40j| can be 
used for the rescaled master equation. Accordingly, we 
look for the probability P{n, t) = P{Nq, t) in the form of 
the WKB ansatz 


P{q,t)=exp[-NS{q)] (4.6) 

where S{q) is a deterministic state function known as the 
action. Intuitively, this approximation expresses the as¬ 
sumption that the probability of occurrence of extreme 
events, such as extinction, lies in the tail of the PDF, 
which falls away steeply from the steady state. Substitut¬ 
ing (14.6|) into the rescaled master equation (12.41) . which 
contains terms of the form Wr{q — r/N), and Taylor- 
expanding terms such as S{q — r/N) around q, we obtain 
to leading order a Hamilton-Jacobi equation H{p, q) = 0 
( 4 ^ . with Hamiltonian 


Hip, 9 ) = X! ^’■('3') [exp(rp) - 1] 


R 

= i?o 9 [exp(ap)-l]-h^g^[exp(-dp)-l]-bg[exp(-p)-l]. 


(4.7) 


Here q is the coordinate, and p = S'{q) is the conjugate 
momentum. The mean-field dynamics can be found by 
writing the Hamilton’s equation q = dpH along the path 
p = 0. This yields the logistic equation as the mean-field 
description of the system dsi), 


djq) 

dt 

= {q) 



(4.8) 


Equation (14.81) has an nontrivial attracting steady state 
at 


if 


q* = -^{a- 1/Ro), 


(4.9) 


aRo > 1 . 


(4.10) 


Note that a bifurcation occurs at i?o = l/o- This implies 
that the population can maintain a long-lived metastable 
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state for a > 1, even if i?o < 1- Going back to the 
mean number of individuals n, the logistic mean-field rate 
equation (14.8|) reads 


dn 

dt 



where 


_ p 1 oA - 7 

r = a/tQ — 1 =- 

d 


and 


K = Nq, 


2{aX — 7 ) 
dfj, 


(4.11) 


(4.12) 


(4.13) 


are the intrinsic growth rate and the carrying capacity, 
respectively. The mean-field logistic equation is com¬ 
pletely specified by the two parameters r and K. On the 
other hand, the stochastic dynamics is characterized by 
up to five parameters, namely a, d, A, /r and 7 . Relations 
(I4.121) and (14.131) provide a constraint for two of them, 
and the other three are free to take different values while 
keeping exactly the same mean-field logistic equation. 

In order to find the MET, we need to find the opti¬ 
mal path to extinction, which is defined as the nontrivial 
heteroclinic orbit that solves the equation H{q,p) = 0 in 
the phase space (q^p) and connects the state (( 7 *,p = 0 ) 
to the extinction state {q = 0,p = pj). Here p/ is the 
solution to the equation qa{pf ) = 0 , namely, the value of 
the momentum along the optimal path to extinction at 
the point where q vanishes. For the system glD we find 
the optimal path to extinction (activation trajectory) 


daip) = 2 


Rq [exp(ap) — 1] — 1 -I- exp(—p) 
Ro [1 - exp(-dp)] 


(4.14) 


and p/ is the solution of the transcendental equation 


i?o [exp(ap/) - 1] - 1 -I- exp{-pf) = 0. (4.15) 

According to [ 2 ^ , the MET is given by 

T = \/^^^^^®xp(iVAS') exp(A^), (4.16) 

where, taking into account (14.141) . 


1 . h'gjP = ^ 
N 


dRo IRq(i(^(i -(- d) -(- 1 — d 

2{aRo - 1) V N ■ 

(4.17) 


The quantities AS and A(j) can be calculated as follows. 
AS” is the action increment along the extinction path, 
which gives us the logarithm of the mean time to extinc¬ 
tion [2^. Since p = dS/dq, 


AS = 5(0) - S{q,) 


Pa{q)dq = / qa{p)dp. (4.18) 


'PS 


Making use of (I4.14p . we obtain from (14.18^ : 


AS = 5(0) - 5(<7*) 


= 2 


- (1 -I- i?o ^)^‘^ + -^0 


/e”/ 


z[z<^ - 1 ) 


-dz. (4.19) 


For d= 1, this equation yields 


AS = S(0) - S(,.) = ^ + 2 V i - 2 V 


For d = 2, (|4.19l) yields for even a. 


(4.20) 


AS = S(0) - S(,.) = 2 (1 + T) In (^i±r|E(w) 

+ («1) 


1 = 1 


2j 


and for odd a. 


A5 = 5(0) - 5(g,) = 2 ( 1 + — ) In 


1 -I- exi 


P(P/) \ 


{a+l)/2 

2 E 

1=1 


1 - exp[( 2 j - l)p/] 
2j - 1 


(4.22) 


In order to go beyond leading-order calculations, we 
determine A(p = 4>{q = 0) — (j){q = q^,), using its definition 
given in [ 2 ^, 



Hpq{qa,p) + \[q'a{p)] ^Hpp{qa,p) + ^qa{p) [exp(-dp) - 1 ] 
Hp{qa,p) 



(4.23) 


where q^ip) = dqa/dp, and the subscripts on H indicate partial derivatives. Making use of (14.7p and (14.141) . we obtain 
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from (14.231) for d = 1, 


2 2 


1 -f- G 


■In 


aexp[(a + l)pf] — (1 + a) exp(ap/) + 1 


and for d = 2, 


A^ = -^ + iln 
^ 2 2 


a [expipf) - 1] 

aRoexp[pf{a + 3)] - Ro{a + 2) exp[p/(a + 1)] - exp(2p/) + 2(i?o + 1) exp(p/) - 1 
{a'^Ro + 2ai?o — 1) [exp(2p/) — 1]^ 

I- 


(4.24) 


(4.25) 


The formula for the coefficient Ai is given by, see 
Eq. (39) of [H, 


Ai — 




where li are the roots of the equation 

n;;(0)r+i-[l + <(0)]/ + l = 0. 


(4.26) 


(4.27) 


Note that the value of does not depend on d. The 
origin of this behavior lies in the fact that reactions with¬ 
out a linear term in n, such as (I4.1bl) . do not play a role 
in the recursive solution of the master equation for small 
values of n [2^. 

Exact analytic expressions for the MET can be ob¬ 
tained for some specific cases. For example, if a = d = 1 
(birth-competition-death), the reactions are given by 


Here we have used Eq. (31) of [^. It can be shown that 
one root of this equation is always ^ = 1. We denote this 
root by Iq. Using the fact that Wq(0) = Rq, and dividing 

by ^ — 1, we need to solve the equation ^i?o(l + ^ H-+ 

Z“-1) -1 = 0, i.e.. 


1 + 


+ U = ^. 

Uq 


(4.28) 


X A 2X, 
2X A X, 

X A0. 


Equation (14.161) yields for the MET, 


(4.32a) 

(4.32b) 

(4.32c) 


For 0 = 1, (14.281) has a single root, and therefore (14.261) 
simplifies to = l/(i?o — !)• For o = 2, (14.281) has 

two roots, and we find 


1 / TT 

r = -1 




, 3/2 


■ exp 


2/V 1 - 


1 -I- In i?o 


= 


3v'i?2 + 4i?o--Ro-4' 


(4.29) 


7 V fv (i?o - ly r*' V* Ro 

recovering the result given by Eq. (70) in 


(4.33) 


For 0 = 3, the polynomial of (14.281) is of third order, and 
its solution yields 


For 0 = 2 and d = 1, the reactions are again of birth- 
competition-death type. 


.(a= 3 ) (c -1- 2)(4 - c){y -2c + +4c+ 16) 

X ^ 3X, 

(4.34a) 

' ^ 3(c2 - 8c - 8)(c4 - 8c2 -h 64) 

(4.30) 

2X A X, 

(4.34b) 

where 

X A0. 

(4.34c) 


c = 


3VS^/3I^TURyT^ + 7Ro + 27 


i?n 


1/3 


(4.31) 


In this case, (j4.16l) yields for the MET, 

^3-\/i?o(d?o + 4) -l- Ro + 4^ Rq -I- 4 -|- ^/Ro{Ro 


1 


j \ N 


2(2i?o- l)2(i?o + 4) 


exp(2A^AS'), 


(4.35) 


where 


AS' = 


■ In 


Rq \ + a/ Ro{Ro + 4) 


{3Ro - 2)^Ro{Ro + 4) + 3i7g + 4i7o - 2 
^i?0 + \/ RoiRo + 4)^ 


(4.36) 
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The difference between results (14.331) and (14.351) is due the value of a. Figure |8] shows that increasing a by one unit 
increases the MET by several orders of magnitude. One can show that T(a = 2)/T(a = 1) ^ exp(A^) as Rq tends to 
infinity. 

Finally, we consider the case a = 1 and d = 2, 

X ^ 2X, 

2X A 0, 

X A 0. 

In this case the reactions are of birth, annihilation, and death type. From (I4.16P we find 

_ 2 [T 

7ViV(A-l)2(i?o + l)i/2 

recovering the result already obtained in [2^. In Fig.[8]we plot the cases o = d = 1 and a = 1, d = 2. The fact that for 
d = 2 we have annihilation, rather than competition as for d = 1, reduces the MET as one would expect. If we take 
i?o to infinity and compare the cases of d = 2 and d = I, both with a = 1, we find r(d = 2)/r(d = I) ~ exp(—2A^ln2). 
For both panels in Fig. |5] we observe that the MET increases very fast with the basic reproductive number Rq. The 
comparison between simulations and analytic results are in general good, except when Rq tends to the critical value 
given in (14.101) . Indeed, the WKB theory breaks down if the barrier IS.S tends to zero. This happens if p/ tends to 0. 
Considering (I4.15|) . the limit p/ —)> 0 implies Rq —)> 1/a. 


exp / 2N 


l + -)ln 
ito 


l±Jk) +1 - J- 

2Ro J Rq 


(4.38) 


(4.37a) 

(4.37b) 

(4.37c) 


V. CONCLUSIONS 

We have adopted an individual-based formulation to 
describe the random dynamics of finite-sized populations. 
Specifically, we have analyzed in detail various possible 
microscopic scenarios that all give rise to the same macro¬ 
scopic population-level model, namely the Verhulst or lo¬ 
gistic population growth equation. We have shown that 
for birth and competition interactions between individu¬ 
als, X A (a -I- 1)X, 2X A X, the population does not 
become extinct, regardless of the value of the parameters. 

If competition leads to annihilation of the competitors, 
xA (a-l-l)X, 2XA0, the ultimate fate of the popula¬ 
tion depends on whether the kinetics is parity conserving 
or not. The parity of the total number of particles is pre¬ 
served in the even-offspring case. This implies that the 
population persists if a is even and uq is odd, because 
the absorbing state is inaccessible. On the other hand, if 
a and uq are both even or if a is odd, the absorbing state 
is accessible and the population becomes extinct. It is 
worth noting that these kinetic rules can be implemented 
as dynamical lattice models or interacting particle sys¬ 
tems, for example as a contact process or a branching- 
annihilating random walk (BARW) (49l - l55j| and that par¬ 
ity conservation, or the lack thereof, also plays a crucial 
role in the dynamics of these spatially extended systems. 
They can display a nonequilibrium transition from a non¬ 
trivial fluctuating steady state to an absorbing state with 
no fluctuations [^. This transition belongs to different 
universality classes for parity-conserving and nonparity- 
conserving models [5^. The most prominent member of 
the first class is the BARW with an even number of off¬ 
springs. The dynamics of BARWs with even and odd 


number of offsprings have been analyzed in detail in . 

For those cases where the population persists, we have 
obtained analytic expressions for the generating function 
and the PDF in the stationary state. In particular, we 
have determined the mean of the PDF and its coefficient 
of variation. For those cases where the population be¬ 
comes extinct, we have calculated the MET and have 
explored its dependence on the microscopic parameters. 
All our analytical results have been compared with nu¬ 
merical simulations, showing good agreement. 

Our results provide further evidence for the advantages 
of individual-based models. They demonstrate that the 
microscopic details of random events at the level of the in¬ 
dividuals lead to differences in the behavior of the system 
at the population level. In the case that the population 
persists, the characteristics of the stationary PDF depend 
on the features of the microscopic model. To illustrate 
this fact, we have focused on the coefficient of variation. 
Our results show, see Figs. [3] and [SJ that an increase in 
a, the number of offsprings, and d, the number of indi¬ 
viduals removed due to competition, leads to an increase 
in the variability of the population for the same value 
of the macroscopic parameter, the carrying capacity JV. 
Measuring the coefficient of variation of a population for 
a given value of JV provides therefore a means of drawing 
inferences about the microscopic details of the birth and 
competition processes. 

Similarly, in the case that the population becomes ex¬ 
tinct, the MET depends sensitively on the microscopic 
details of the model, as illustrated by Figs. HI [71 and HI 

For example, for the model X A (a -I- 1) X, 2 X A 0, 
we find that if no and a are even, then the MET be¬ 
comes significantly larger for the same carrying capacity 
as a increases, see Fig. HI In contrast, the MET becomes 
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Ro 

FIG. 8 : Mean time to extinction for birth-and-death reac¬ 
tions. In both cases we set /i = 0.1, A = 10, d = 1, = 100 

and modify 7 to vary Rq. In panel a) we compare different 
values for the MET for different d and the same a. In panel b) 
we compare different values for the MET for different a and 
the same b. Simulations (symbols) have been performed up 
to time 10 ®, and mean values are obtained by averaging over 
4 X 10^ realizations. Solid curves correspond to exact analytic 
results given by (14.331) . (14.3511 and (14.381) 


significantly smaller for the same carrying capacity as a 
increases if a is odd, see Fig. 0 Extinction is always 
the ultimate fate for the birth-competition-death model, 
X A- (a-l-l)X, 2X A (2 —d)X, X A 0. Figure [8] demon¬ 
strates strikingly the sensitive dependence of the MET 
on the microscopic details of the system. Measuring the 
MET for laboratory populations with a given basic re¬ 
productive number provides therefore a means of draw¬ 
ing inferences about the microscopic details of the birth, 
death, and competition processes. Our results also im¬ 
ply that assessing the extinction risks and survival times 
of natural populations requires an understanding of the 
microscopic details of the processes that occur in the sys¬ 
tem and should not be based solely on phenomenological 
models. 

There are many other possibilities that lead to the lo¬ 
gistic equation, for example if we consider two different 

birth reactions simultaneously, such as, X A 2 X and 
X A 3X. Another possibility consists in considering 
schemes with four or even more reactions. All these situ¬ 
ations can be analyzed in the same manner and with the 
same techniques as used here. A further intriguing pos¬ 
sibility that deserves study are reactions schemes where 
the number of offsprings, a, and the number of individu¬ 
als eliminated by exclusive competition, d, fluctuate ran¬ 
domly between several values. 
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